---
title: "Power Analysis"
date: "8/19/2019"
output: pdf_document
header-includes:
 \usepackage{float}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE, fig.pos = 'H')
library(haven)
library(dplyr)
library(DeclareDesign)
library(ggplot2)
```

First, using the existing data, we extract the standardized effect sizes for the difference-in-means, and the interaction effect.

```{r, messages = FALSE, warning = FALSE}
library(haven)

pipelines_original <- read_dta("pipelines_fullvar.dta")
pipelines <- pipelines_original[,c('accessed', 'disagreement')]

simple_ate <- mean(subset(pipelines_original, Z == 1)$accessed) -
  mean(subset(pipelines_original, Z == 0)$accessed)

pooled_sd <- sqrt((sd(subset(pipelines_original, Z == 1)$accessed) +
  sd(subset(pipelines_original, Z == 0)$accessed))/2)
effect_size <- simple_ate/pooled_sd

# standardized interaction effect
interaction_effect <- as.numeric(coef(lm(accessed ~ Z*disagreement,
                      data=pipelines_original))[4])
interaction_effect_size <- interaction_effect/(sqrt(sd(pipelines_original$accessed)))
```

The design declaration for the first model is fairly straightforward: a difference-in-means using the observed ATE and noise that is distributed normally with $mu = 0$ and $sigma^2$ = 1. A random number seed is used for reproducibility.

```{r}
# set seed to preference
set.seed(08192019)

library(DeclareDesign)

design_model1 <- declare_population(N = N, noise = rnorm(N)) +
  declare_potential_outcomes(Y ~ effect_size * Z + noise) +
  declare_estimand(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_assignment(prob = 0.5) +
  declare_reveal(Y, Z) +
  declare_estimator(Y ~ Z, model = lm_robust, term = "Z", estimand = "ATE"
)
```

With the second design declaration, for diagnosing the design at larger sample sizes than the existing observed sample size, it is necessary to assume a distribution for the included pre-registered covariate, `disagreement`, which represents constituency disagreement. In this case, it is normally distributed using the observed mean in the dataset and its observed standard deviation, or $N \sim (-0.161, 0.265)$. This is a plausible assumption given the existing distribution of the covariate, as shown in Figures 1 and 2, though the results are substantively similar using a different distribution. 

\newpage

```{r, echo = FALSE, fig.cap="Distribution of observed constituency disagreement", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
library(ggplot2)

ggplot(pipelines, aes(disagreement)) +
  geom_histogram(binwidth = 0.2, color="black", fill="white") +
  scale_x_continuous(limits=c(-2,2)) +
  stat_function(fun = function(x) dnorm(x, mean = mean(pipelines$disagreement), sd = sd(pipelines$disagreement)) * 331 * 0.2, color = "darkred", size = 1) +
  theme_bw() +
  labs(x="Constituency disagreement scale", y = "Count")
```

```{r, echo = FALSE, fig.cap="Q-Q plot of constituency disagreement", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
qplot(sample = disagreement, data = pipelines) +
  theme_bw() +
  labs(y = "Observed", x = "Theoretical")
```

\newpage

The design declaration for the second model uses the same observed ATE and parameters for the noise as the first model, and then uses the observed standardized interaction effect size. The estimand is a difference in means acounting for the continous interaction effect.

```{r}
design_model2 <- declare_population(N = N,
                  disagreement = rnorm(N, mean(pipelines$disagreement), sd(pipelines$disagreement)),
                  noise = rnorm(N)) +
  declare_potential_outcomes(Y ~ effect_size * Z + effect_size * disagreement +
                               interaction_effect_size * disagreement * Z + noise) +
  declare_estimand(ATE = mean(Y_Z_1 - Y_Z_0), interaction = interaction_effect_size) +
  declare_assignment(prob = 0.5) +
  declare_reveal(Y, Z) +
  declare_estimator(
    Y ~ Z + disagreement + Z * disagreement,
    model = lm_robust,
    term = c("Z", "Z:disagreement"),
    estimand = c("ATE", "interaction")
  )
```

Each of the designs is assesed at $N = 331, 500, 1000, 3000$ for 1,000 simulations.

```{r}
designs_model1 <- redesign(design_model1, N = c(331, 500, 1000, 3000, 5000))
designs_model2 <- redesign(design_model2, N = c(331, 500, 1000, 3000, 5000))
simulations_model1 <- simulate_design(designs_model1, sims = 1000)
simulations_model2 <- simulate_design(designs_model2, sims = 1000)
```

Figure 3 shows the power analyses at each of the designated sample sizes for the difference-in-means estimand.

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.cap = "Power Analysis for Model 1", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
library(dplyr)

summary_df_1 <-
  simulations_model1 %>%
  mutate(p.value_one_sided =
           ifelse(estimate < 0, p.value/2, (1-p.value)/2)) %>%
  group_by(N, estimand_label) %>%
  summarize(power = mean(p.value_one_sided <= 0.05),
            mean_estimate = mean(estimate))

ggplot(summary_df_1, aes(N, power, group = estimand_label, color = estimand_label)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0,1)) +
  labs(colour = "Estimand", y = "Power", x = "Sample Size")
```

\newpage

Figure 4 shows the power analyses at each of the designated sample sizes for both estimands in the second model.

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.cap = "Power Analysis for Model 2", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
summary_df_2 <-
  simulations_model2 %>%
  mutate(p.value_one_sided =
           ifelse(estimand_label == "ATE" & estimate < 0, p.value/2, ifelse(estimand_label == "ATE" & estimate > 0, (1-p.value)/2, ifelse(estimand_label == "interaction" & estimate > 0, p.value/2, (1-p.value)/2)))) %>%
  group_by(N, estimand_label) %>%
  summarize(power = mean(p.value_one_sided <= 0.05),
            mean_estimate = mean(estimate))

ggplot(summary_df_2, aes(N, power, group = estimand_label, color = estimand_label)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0,1)) +
  labs(colour = "Estimand", y = "Power", x = "Sample Size")
```

Using a simulation-based approach, it is clear that the estimates of both the ATE and the interaction effect are small enough that even if the sample size were much larger, it is unlikely for the result to have been statistically significant.

An alternative approach would be to vary the effect size and hold the sample size constant at 331 to see, given the observed sample size, how large the ATE and interaction effect would have to have been to have estimated a statistically singificant result. This strategy requires no additional assumptions about the constituency disagreement variable.

```{r}
design_model1 <- declare_population(N = 331, noise = rnorm(331)) +
  declare_potential_outcomes(Y ~ effect_size * Z + noise) +
  declare_estimand(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_assignment(prob = 0.5) +
  declare_reveal(Y, Z) +
  declare_estimator(Y ~ Z, model = lm_robust, term = "Z", estimand = "ATE"
)

design_model2 <- declare_population(N = 331,
                  disagreement = pipelines$disagreement,
                  noise = rnorm(331)) +
  declare_potential_outcomes(Y ~ effect_size * Z + effect_size * disagreement +
                               interaction_effect_size * disagreement * Z + noise) +
  declare_estimand(ATE = mean(Y_Z_1 - Y_Z_0), interaction = interaction_effect_size) +
  declare_assignment(prob = 0.5) +
  declare_reveal(Y, Z) +
  declare_estimator(
    Y ~ Z + disagreement + Z * disagreement,
    model = lm_robust,
    term = c("Z", "Z:disagreement"),
    estimand = c("ATE", "interaction")
)
```

In this case, we vary the effect size starting with the observed effect sizes and then according to the typical values for Cohen's D, where 0.2 is considered to be a small effect, 0.5 is considered to be a medium-sized effect, and 0.8 is considered to be a large effect.

```{r}
observed_effect_size <- effect_size
observed_interaction_effect_size <- interaction_effect_size
rm(effect_size)
rm(interaction_effect_size)

designs_model1 <- redesign(design_model1,
                           effect_size = c(observed_effect_size, -0.2, -0.5, -0.8))
designs_model2 <- redesign(design_model2,
                           effect_size = c(observed_effect_size, -0.2, -0.5, -0.8),
                           interaction_effect_size = c(observed_interaction_effect_size, 0.2, 0.5, 0.8))
simulations_model1 <- simulate_design(designs_model1, sims = 1000)
simulations_model2 <- simulate_design(designs_model2, sims = 1000)
```

Figures 5-9 show the results. The study is sufficiently powered to detect small, medium, and large effect sizes.

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.cap = "Power Analysis for Model 1", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}

summary_df_1 <-
  simulations_model1 %>%
  mutate(p.value_one_sided =
           ifelse(estimate < 0, p.value/2, (1-p.value)/2)) %>%
  group_by(effect_size, estimand_label) %>%
  summarize(power = mean(p.value_one_sided <= 0.05),
            mean_estimate = mean(estimate))

summary_df_1$effect_size <- case_when(
  summary_df_1$effect_size == "observed_effect_size" ~ "Observed Effect Size",
  summary_df_1$effect_size == "-0.2" ~ "Small",
  summary_df_1$effect_size == "-0.5" ~ "Medium",
  summary_df_1$effect_size == "-0.8" ~ "Large"
)

summary_df_1$effect_size <- ordered(summary_df_1$effect_size, levels = c("Observed Effect Size", "Small", "Medium", "Large"))

ggplot(summary_df_1, aes(effect_size, power, group = estimand_label, color = estimand_label)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0,1)) +
  labs(colour = "Estimand", y = "Power", x = "Effect Size")
```

```{r, echo = FALSE, message = FALSE, warning = FALSE}
summary_df_2 <-
  simulations_model2 %>%
  mutate(p.value_one_sided =
           ifelse(estimand_label == "ATE" & estimate < 0, p.value/2, ifelse(estimand_label == "ATE" & estimate > 0, (1-p.value)/2, ifelse(estimand_label == "interaction" & estimate > 0, p.value/2, (1-p.value)/2)))) %>%
  group_by(effect_size, interaction_effect_size, estimand_label) %>%
  summarize(power = mean(p.value_one_sided <= 0.05),
            mean_estimate = mean(estimate))

summary_df_2$effect_size <- case_when(
  summary_df_2$effect_size == "observed_effect_size" ~ "Observed",
  summary_df_2$effect_size == "-0.2" ~ "Small",
  summary_df_2$effect_size == "-0.5" ~ "Medium",
  summary_df_2$effect_size == "-0.8" ~ "Large"
)

summary_df_2$effect_size <- ordered(summary_df_2$effect_size, levels = c("Observed", "Small", "Medium", "Large"))

summary_df_2$interaction_effect_size <- case_when(
  summary_df_2$interaction_effect_size == "observed_interaction_effect_size" ~ "Observed",
  summary_df_2$interaction_effect_size == "0.2" ~ "Small",
  summary_df_2$interaction_effect_size == "0.5" ~ "Medium",
  summary_df_2$interaction_effect_size == "0.8" ~ "Large"
)

summary_df_2$interaction_effect_size <- ordered(summary_df_2$interaction_effect_size, levels = c("Observed", "Small", "Medium", "Large"))

interaction_subset <- subset(summary_df_2, estimand_label == "interaction")
ATE_subset <- subset(summary_df_2, estimand_label == "ATE")
```

\newpage

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.cap = "Power Analysis for ATE Estimand for Model 2, Varying Both Effect Sizes", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
ggplot(ATE_subset, aes(effect_size, power, group = interaction_effect_size, color = interaction_effect_size)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0,1)) +
  labs(colour = "Interaction Effect Size", y = "Power", x = "Effect Size")
```

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.cap = "Power Analysis for ATE Estimand for Model 2, Varying Both Effect Sizes", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
ggplot(ATE_subset, aes(interaction_effect_size, power, group = effect_size, color = effect_size)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0,1)) +
  labs(colour = "Effect Size", y = "Power", x = "Interaction Effect Size")
```

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.cap = "Power Analysis for Interaction Estimand for Model 2, Varying Both Effect Sizes", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
ggplot(interaction_subset, aes(effect_size, power, group = interaction_effect_size, color = interaction_effect_size)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0,1)) +
  labs(colour = "Interaction Effect Size", y = "Power", x = "Effect Size")
```

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.cap = "Power Analysis for Interaction Estimand for Model 2, Varying Both Effect Sizes", fig.height = 3, fig.width = 5, fig.show = "hold", fig.align= "center"}
ggplot(interaction_subset, aes(interaction_effect_size, power, group = effect_size, color = effect_size)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_y_continuous(limits = c(0,1)) +
  labs(colour = "Effect Size", y = "Power", x = "Interaction Effect Size")
```